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DETERMINING AMPLITUDE LIMITS FOR VIBRATION SPECTRA 

BACKGROUND 

[0001] The present invention relates generally to the field of detecting 
mechanical faults in HVAC equipment. More specifically, the present invention relates 
to a methodology for calculating vibration amplitude limits to detect mechanical faults in 
HVAC equipment such as chillers. 

[0002] Timely detection, diagnosis and repair of mechanical problems in 
machinery such as heating, ventilating and air-conditioning (HVAC) systems is important 
for efficient operation. Chillers are important components of HVAC systems because 
they consume a large fraction of energy in a building and require a large capital 
investment. Severe mechanical faults in chillers typically results in expensive repairs and 
disruptions to the HVAC system during the repair period. Accordingly, chillers are 
generally monitored routinely to detect developing mechanical faults. 

[0003] A common method for detecting and diagnosing mechanical faults is 
vibration analysis. By analyzing vibration data from different positions on a chiller, a 
vibration analyst can detect and diagnose mechanical faults in a machine. Vibration data 
is commonly available as a spectrum. The vibration analyst can determine a machine's 
condition by analyzing the vibration amplitude at different frequencies in the spectrum. 
A vibration analyst typically detects a mechanical fault that requires corrective action 
when the amplitudes in the vibration spectrum exceed "acceptable" limits. Under many 
current approaches, the acceptable limits are specified by rules-of-thumb or from a 
vibration analyst's experience. However, these approaches can be unreliable. For 
example, limits determined from an individual's experience can be incorrect or 
inconsistent. Similarly, limits specified by rules-of-thumb are typically generalized to 
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apply to a large number of chillers and are therefore not likely relevant for some 
particular types of chillers. 

[0004] Accordingly, there exists a need for a method of more accurately 
detecting mechanical faults without having to rely on an individual's knowledge of a 
system or general rule-of-thumb limits. In particular, it is desirable to be able to derive 
the limits from historical data using advanced statistical methods. By using statistics and 
historical data, the estimated limits may be based entirely on the vibration spectrum of 
each type of chiller. This approach easily allows updating of the amplitude limits when 
new vibration data from chillers is collected. In addition, because this approach uses 
statistical methods and not expert knowledge or "rules-of-thumb," it results in more 
consistent limits. Many vendors also use rudimentary statistics such as calculating the 
average and standard deviation of the data to estimate limits. Unfortunately, this 
approach can result in erroneous limits because vibration data usually don't follow a bell 
shaped (or Gaussian) probability distribution. 

[0005] It would be advantageous to provide a method or the like of a type 
disclosed in the present application that provides any one or more of these or other 
advantageous features. The present invention further relates to various features and 
combinations of features shown and described in the disclosed embodiments. Other ways 
in which the objects and features of the disclosed embodiments are accomplished will be 
described in the following specification or will become apparent to those skilled in the art 
after they have read this specification. Such other ways are deemed to fall within the 
scope of the disclosed embodiments if they fall within the scope of the claims which 
follow. 

SUMMARY 

[0006] One embodiment of the present invention relates to a method for 
determining vibration amplitude limits to detect faults in mechanical equipment. The 
method comprises estimating a data probability distribution based on data for the 
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mechanical equipment and utilizing the data probability distribution to calculate the 
vibration amplitude limits. 

[0007] Another embodiment of the present invention relates to a method for 
detecting faults in a chiller based on vibration amplitude limits. The method comprises 
calculating vibration amplitude limits of the chiller using statistics and historical data for 
the chiller, estimating an at least two-dimensional density estimate, and weighting the 
historical data based on when the historical data was generated. The vibration amplitude 
limits are calculated as a function of frequency for an entire frequency spectrum. 

[0008] Another embodiment of the present invention relates to a method for 
determining vibration amplitude limits of a mechanical device. The method comprises 
identifying a mechanical device and a frequency range for a spectrum to be analyzed, 
retrieving vibration spectra comprising individual spectrum for the mechanical device 
and the frequency range, calculating frequency for the individual spectrum, and 
identifying the individual spectrum with the smallest number of frequency lines. In 
addition, the method comprises calculating noise bandwidths and a largest noise 
bandwidth, removing outlier data, calculating conditional kernel density, and calculating 
vibration amplitude limits to detect faults in the mechanical device. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0009] FIGURE 1 is a block diagram illustrating a method of estimating 
vibration amplitude limits according to an exemplary embodiment. 

[0010] FIGURE 2 is a diagram representing a function for removing old data 
from a dataset according to an exemplary embodiment. 

[001 1] FIGURE 3 is a diagram representing frequency interpolation according 
to an exemplary embodiment. 

[0012] FIGURE 4 is a block diagram illustrating a method of detecting outlier 
data from a dataset according to an exemplary embodiment. 

[0013] FIGURE 5 is a diagram representing the application of multivariate 
outlier removal from vibration data according to an exemplary embodiment. 
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[0014] FIGURE 6 is a diagram illustrating a variety of kernel functions 
according to an exemplary embodiment. 

[0015] FIGURE 7 is a diagram comparing kernel density estimates obtained 
using different bandwidth selection methods according to an exemplary embodiment. 

[0016] FIGURE 8 is a diagram illustrating Epanechnikov boundary kernels for 
different values according to an exemplary embodiment. 

[0017] FIGURE 9 is a diagram illustrating the effect of boundary bias on kernel 
density estimation according to an exemplary embodiment. 

[0018] FIGURE 10 is a diagram illustrating the effect of weighting on the alert, 
alarm, and danger confidence levels with diagnostic frequencies for a chiller according to 
an exemplary embodiment. 

[0019] FIGURE 1 1 is a diagram illustrating amplitude limit envelopes for the 
motor- vertical position of a chiller according to an exemplary embodiment. 

[0020] FIGURE 12 is a diagram illustrating amplitude limit envelopes for the 
compressor-vertical position of a chiller according to an exemplary embodiment. 

DETAILED DESCRIPTION 

[0021] Before explaining a number of exemplary embodiments of the invention 
in detail, it is to be understood that the invention is not limited to the details or 
methodology set forth in the following description or illustrated in the drawings. The 
invention is capable of other embodiments or being practiced or carried out in various 
ways. It is also to be understood that the phraseology and terminology employed herein 
is for the purpose of description and should not be regarded as limiting. 

[0022] In general, the method described in this disclosure uses observed data to 
identify the probability distribution for the data and then uses the probability distribution 
to calculate limits for the data. For example, if a user finds that data follows a Gamma 
probability distribution, then the user can fit the data to a Gamma distribution, and then 
calculate the limits from the probability distribution. In many instances, vibration data do 
not follow a single type of distribution, and therefore there is a need to automatically 
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generate the probability distribution from the data itself. A preferred embodiment of the 
invention calculates the probability distribution from the data using the kernel density 
estimation method, and then uses the estimated probability distribution to calculate the 
vibration limits. 

[0023] An exemplary methodology for calculating vibration limits from 
historical data will now be discussed. This methodology is shown in FIGURE 1. Step 1 
requires user input while the rest of the steps are executed automatically. 

[0024] At step 1, a vibration analyst specifies the chiller model for which the 
limits will be calculated as well as the frequency range of the spectrum. The frequency 
range of the spectrum is from Frange, i ow to F range , h i Hz. The analyst also specifies the 
electrical line frequency (e.g., 60 Hz in North America, 50 Hz in Europe) and the 
importance levels for different frequencies. The importance levels for different 
frequencies is described in greater detail below. 

[0025] At step 2, a historical database is queried to retrieve all spectra for the 
specified chiller make and model and the selected frequency range. The set denotes 
the collection of n spectra retrieved from the historical data and the /* member of S is 
denoted as Sj. Each spectrum Sj consists of frequency and amplitude data. 

[0026] At step 3, a motor frequency M/j (/ = 1, ...,«) is estimated for every 
spectrum Sj and the frequency axis of Sj is scaled by its motor frequency, Af/j. Scaling 
the frequency axis by the motor frequency reduces the variation in the spectra due to 
differences in motor speeds. The motor frequency for every spectrum Sj is calculated by 
first searching for the highest amplitude in the range 57-60 Hz (60 Hz is the electric line 
frequency in North America) and then using interpolation or another method to refine the 
frequency location of the peak. For example, an interpolation method may be used such 
as described by Technical Associates, Concentrated Vibration Signature 
Analysis and Related Condition Monitoring Techniques (2002). FIGURE 2 
illustrates this method. According to an exemplary embodiment illustrated in FIGURE 3, 
while locating the motor frequency from the vibration spectrum,^, may be the frequency 
value of the highest peak in the range of 57-60 Hz. lfA p is the amplitude at f p9 f s and A s 
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are the frequency and amplitude of the line having an amplitude closest to A P9 and Af is 
the spacing between adjacent frequency lines, then the interpolated frequency ft may be 
calculated as: 

yA s +Apj 

J -\f P -Afciffi<f P 
[0027] At step 4, the spectrum having the smallest number of frequency lines in 
S is found. If the number of frequency lines in spectrum Sj is m h then m = min (mj) . The 

l£/£fi V 7 

value of m results in the calculation of the largest noise bandwidth Ntw in S. 

[0028] At step 5, the largest noise bandwidth Nbw in S is calculated as: 

Frange, hi - Frange, lo 



Mw = 1.5x 



mxMf 



(2) 



where M/ is the average value of the machine frequency for n spectra: 

M/^YM/j (3) 

According to an exemplary embodiment, the factor 1.5 is utilized assuming that a 
Hanning window was used for calculation of the spectrum from its time waveform. 

[0029] At step 6, a loop begins for calculating the vibration amplitude limits at 
different frequencies. Although the limits can be calculated at any frequency, according 
to a preferred embodiment, the limits are calculated at frequencies spaced less than or 
equal to Ntw /3 to have a good resolution for the amplitude limit spectrum. Using results 
from step 4, there are 2m+l frequency values where the limits are calculated. Thus, step 6 
implements a loop counter, i, from 0 to 2m. 

[0030] At step 7, a frequency value at which the amplitude limits will be 
estimated is calculated: 

_ Frange, lo + iX Nbw/ 3 

f " M, (4) 
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where fi is scaled by the average motor frequency. 

[0031] At step 8, after the frequency value has been specified at step 7, 
vibration data are collected in the frequency range fi - 2Nbw to fi + 2Ntw from each 
spectrum Sj. In general, step 8 involves collecting vibration data in the range 
f±(cxNbw) , where c is an integer >2. The preferred value c = 2 ensures that all data 
within the noise bandwidth is included in the density calculation and the computational 
load is at a minimum. A value of c < 2 would result in elimination of data just outside 
the noise bandwidth. Recall that step 3 had scaled the frequency axis of each spectrum 
by its motor frequency. The data collected in step 8 has two dimensions: frequency and 
amplitude. 

[0032] At step 9, outlier detection and removal for the two dimensional data 
collected from step 8 is implemented. Outlier removal is helpful to reduce the effect of 
unusual observations in limits calculation. Outliers are observations that appear to be 
inconsistent with the majority of data in a dataset. These observations are preferably 
removed from the data because they may corrupt the data analysis and can produce 
misleading results. For example, the average for the dataset, X= {10, 9, 1 1, 8, 12, 50, 1 1, 
30, 8}, is 16.6, even though the majority of the data are centered around 10. The 
observations 30 and 50 push the average to a higher value that suggests that most of the 
observations are centered around 16, which could be a misleading result. These 
observations also appear to be much larger than the majority of the data, and thus are 
called outliers. The average becomes 9.86 when these extreme observations are removed 
from the dataset, a result that is consistent with the majority of observations. 

[0033] Many types of outlier detection procedures are known in the art. 
According to an exemplary embodiment, outlier removal using kernel density estimation 
may be utilized. According to a preferred embodiment, the sequential Wilk's 
multivariate outlier detection test may be utilized to determine outliers, as explained by 
C. Caroni and P. Prescott, Sequential Application of Wilk 's Multivariate Outlier Test, 
Appl. Statistics 41, 355-64 (1992), hereinafter referred to as the "C and P method". 
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The C and P method is well suited to the present embodiment because the data is two- 
dimensional (frequency and amplitude). The detailed procedure for detecting outliers 
using the C and P method will now be described. 

[0034] The algorithm for detecting multivariate outliers in data is presented in 
FIGURE 4. Given a dataset X= {x\, xi, . . . , x n }, with n observations in d dimensions, 
and an estimate of the upper bound r for the number of outliers, the procedure presented 
in FIGURE 4 identifies possible outliers in the dataset X. The descriptions of the steps of 
the procedure are presented below. At step 21, the number of outliers is initialized to 
zero by setting « out = 0. At step 22, the average x of the observations for dataset X is 
computed as: 



where xj is a member of the dataset X and n equals the number of observations in X. 

[0035] At step 23, the matrix of sums of squares of outer products, A h for the 
dataset X is computed. The matrix sum of squares is determined as: 



[0036] At step 24, a test is performed to determine whether the matrix sum of 
squares is zero (A t = 0). If this matrix is zero, then all of the observations in dataset X 
have the same value and there are no outliers in the remaining observations in dataset X. 
To prevent division by the determinant of a zero matrix at step 25, the process moves to 
step 31 when the matrix sum of squares in step 24 equals zero. Otherwise, the process 
advances to step 25. 

[0037] At step 25, the i th extreme observation x eJ in the dataset Xis found. The 
most extreme observation x e j can be identified where its removal results in the smallest 
value of W/ = where A iA is the matrix obtained by removing the most extreme 

observation in the previous step. The notation \A\ refers to the determinant of the matrix 
A. 




(5) 



4 = X(*/ -*X*/ _ *) r 



(3) 
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[0038] At step 26, the i th extreme Wilk's statistic A is computed. The Wilk's 
statistic is determined as: 



Di = min 



\Ai\ 



Ai- 



(7) 



[0039] At step 27, the / th critical value i is computed. The C and P method 
utilizes the following relation for determining the critical value: 



fl d Z7 

L n-d-i 



+ 1)) 



(8) 



where Fd,n-d- i, a - a /(« - 1 + 1» is the 100x(l-o/(n - i + l)) th percentage point for the F 
distribution with d and {n-d-i) degrees of freedom, a is the probability of detecting 
outliers when there are no outliers, and d is the dimensionality of the data. 

[0040] At step 28, a test is performed to determine if the i th critical value A is 
less than the critical value y. (D. < y i ), determined in step 27. At step 29, the number of 
outliers n out is set equal to i (n out = i). At step 30, extreme observation x e>i is removed 
from the dataset X. After removing the extreme element x e j, the number of observations 
in the dataset X is n - L If i equals r, then the process proceeds to step 3 1 . Otherwise, i is 
incremented and the process returns to step 21. At step 31, the extreme observations that 
are outliers { x e ,u Xe,2, ■ • • , Xe.nout } are identified. The first n out extremes identified at step 
25 are considered outliers. However, not all the extreme values determined at step 25 are 
outliers. 

[0041] The C and P method is designed for Gaussian data. However, vibration 
data does not always follow a Gaussian distribution. Despite this limitation, the C and P 
method may still be used to process the vibration data because the approach enables the 
detection of observations that have very high amplitudes compared to a majority of the 
data. 

[0042] The maximum number of expected outliers in the data, r, is another 
parameter of the C and P method. A large value for r means that more data points are 
expected to be outliers, which results in a large computational effort. According to an 
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exemplary embodiment, r may be set to be the largest integer less than 0.2 x w> where n is 
the number of spectra in the set S to reduce the computational load. 

[0043] The vibration database contains data from acceptable as well as faulty 
machines. Outlier removal is performed to remove the unusual values that are a result of 
severe faults and/or measurement errors. After outlier removal, the database may still 
contain data from machines that have large vibration levels, but these data should be 
included for calculation of amplitude limits in order to observe the full range and 
distribution of vibration amplitudes. The application of the C and P method is presented 
in FIGURE 5 for sample data from a chiller at/= 0.5. 

[0044] Referring back to FIGURE 1, at step 10, the conditional kernel density 
p(x\F = ft) is calculated, where x is the amplitude and F represents the frequency axis. 

An exemplary procedure for calculating p{x\F = fi) from the two dimensional data is 

discussed below. At step 10, the cumulative probability function, P(x\F = fi) , is also 

calculated by numerically integrating the density p(x\F = fi) from zero to jc: 

P (A F = fi) = £ p( z \ f = fi) dz - The amplitude axis is divided into a fine grid of 1000 
points and P(x\F = fi) is calculated numerically at each of these grid points. This 
procedure makes calculating the inverse of P(x\F = f) at step 1 1 relatively easy. 

[0045] At step 1 1, the value x$ is calculated such that P(x\F = f)=0 9 where 0 

is a specified probability cutoff. Thus, xp is the 100 x jS % confidence limit for the 
vibration amplitude. Typical values for (3 are 0.95, 0.99, etc. An exemplary methodology 
for choosing the value for /3 is discussed below. The value of xp can be calculated by a 
look-up table of P(x\F = f) and x grid values from step 10. After the completion of step 

1 1, the process returns back to the beginning of the "For" loop at step 6 and the amplitude 
limits for a new frequency value are calculated. A spectral amplitude limit envelope is 
obtained at the termination of the "For" loop. The vibration analyst can compare the 
spectrum for a new machine with the calculated limit envelope to determine if any 
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amplitude in the new spectrum violates the limits. If any of the amplitudes violate the 
limits, then a machine fault is present. In addition, the frequency location of a limit 
violation provides the vibration analyst with clues about the possible cause of the fault. 

[0046] The kernel density estimation method briefly discussed at step 10 above 
will now be described in greater detail. The kernel density method is described in detail 
by M.P. Wand and M.C. Jones, Kernel Smoothing (1995). The kernel density 
estimation method is also known as the Parzen windows method for density estimation. 
It attempts to estimate an unknown probability density for a given dataset X The 
probability density estimate at a point x for a one-dimensional dataset with n data points 
is given by: 

where, xj is the / h observation of the dataset X, h is called the bandwidth that 
characterizes the spread of the kernel, and k (•) is a kernel density function that is 
symmetric and satisfies the condition: 

^K(u)du = 1 

Some common univariate kernel density functions, k (•), are listed in FIGURE 6. The 
notation 1 <ij means: 

A fl if|w|<l 
fl"M~[0 otherwise (H ' 

[0047] To calculate the kernel density for two-dimensional vibration data (e.g., 
frequency and amplitude directions), an expression should be utilized for the kernel 
density estimation in two or more dimensions. In general, Equation (9) is extended to d 
dimensions by replacing the univariate kernel k (u) by a rf-dimensional kernel K (w), and 
the bandwidth h by a bandwidth matrix H. Then, the rf-dimensional kernel density 
estimate is written as: 

-11- 

001.1581165 



(9) 



(10) 
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p(x) = ±t\H\- U2 K(H- u \x-xj)) (12) 

where | ■ | denotes the matrix determinant. There are two known ways of constructing the 
multivariate density in Equation (12) from univariate kernel functions. One option is to 

use a spherically symmetric multivariate kernel K 5 (w) = c K d K ((m'm) 1/2 ) , where c K d is a 

constant dependent on the univariate kernel k (•), and the dimensionality. A second 

option is to use a multivariate product kernel K p (u) = k(ih) . Although K s (•) and 

K p (0 are usually different, the difference between the kernel densities estimated using 
K s (•) and K p (•) is typically very small. The product kernel K p (•) is easier to implement 
in software, particularly when boundary correction (described below) is not required in 
all dimensions. Thus, the product form of the multivariate kernel is used to calculate the 
density according to a preferred embodiment. 

[0048] For the two-dimensional vibration data, F and x denote the frequency 
and the amplitude dimensions, respectively. Then, the expression for the kernel density 
using the product kernel in these two dimensions is written as: 



nj^hfhx \ hf ) 



fr-r \ 



1 ^ 1 (F-fj\ \ x-x 



K 



h x 



(13) 



In Equation (13), hf and h x are bandwidths of the frequency and amplitude kernels, and fj 
and Xj are the frequency and amplitude values of the / h observation in the dataset, 
respectively. 

[0049] To calculate the vibration amplitude limits at a given frequency value F 
= /, the conditional probability density p(x\F = fi) is first calculated from the joint 

probability density p(x 9 F), Conditional density is then used to calculate the limits. The 
conditional probability density of the vibration amplitude (x) for a given frequency value 
(F - f) is calculated as: 
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p(*\F = f) = - 



(14) 



00 



jp(x,F = f)dx 



[0050] The kernel density estimation method has two important design 
parameters including the choice of the kernel function and its bandwidth. First, although 
the choice of the kernel function is not particularly important based on the efficiencies of 
different kernel functions, the choice is useful from a computational perspective for large 
data sets. For example, a Gaussian kernel requires a much larger computational effort 
than the Epanechnikov kernel. Also, a Gaussian boundary kernel requires the numerical 
calculation of integrals while the expression for the Epanechnikov boundary kernel is 
obtained analytically. 

[0051] The Epanechnikov kernel has the lowest asymptotic mean integrated 
squared error (AMISE) in estimating the kernel density as well as a low computational 
cost. Thus, the Epanechnikov kernel is the preferred kernel for estimating the probability 
density. 

[0052] The second useful parameter is the kernel bandwidth. The kernel 
bandwidth calculation methods can be divided into two broad categories: (i) first 
generation rules-of-thumb, and (ii) second generation methods that produce bandwidths 
close to the optimal value. 

[0053] A rule of thumb (ROT) using the sample variance and the inter-quartile- 
range (IQR) may be used to estimate the bandwidth of the kernels. According to an 
exemplary embodiment, a robust estimate of the kernel width, h, is, 



The quantity a is the sample standard deviation. The inter-quartile-range (IQR) is 
computed as the difference between the 75th percentile (Q3) and the 25th percentile (Ql) 
of the data. Thus, IQR = Q3 - Ql. The quantity represents the width estimate based on 



h =Q.9xAxn 



-1/5 



(15) 



where, 




(16) 
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the inter-quartile range and the factor 1.349 is the inter-quartile-range for a Gaussian 
distribution with zero mean and unit standard deviation. This expression is valid for 
univariate data; 

[0054] Second generation methods are designed to produce bandwidths that are 
close to the optimum but require a larger computational effort. According to an 
exemplary embodiment, the Sheather- Jones bandwidth selection method may be used to 
calculate the kernel bandwidth. This method calculates the bandwidth by using a series of 
"direct-plug-in" equations and is described by S.J. Sheather and M.C. Jones, A 
Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation, J. 
Royal Statistical Soc. Ser. B, v. 53, 683-90 (1991). 

[0055] The Sheather- Jones direct-plug-in (SJ-dpi) method calculates the kernel 
bandwidth by using fourth and higher order derivatives of the kernel function and 
requires them to be non-zero and finite. Because the third and higher-order derivatives of 
the Epanechnikov kernel function presented in this disclosure are zero, the SJ-dpi method 
cannot be used directly to calculate the bandwidth for the Epanechnikov kernel. To solve 
this problem, an exemplary method uses a result called the theory of equivalent kernels, 
which is described by D.W. Scott, Multivariate Density Estimation: Theory, 
Visualization and Practice (1992). By using this theory, the bandwidth calculated 
using one kernel function can easily be transformed to another kernel function. For 
example, if the bandwidth is calculated assuming a Gaussian kernel, then the equivalent 
bandwidth for the Epanechnikov kernel is: 



Because the Gaussian kernel function has infinite number of non-zero and finite 
derivatives, the method calculates the bandwidth using the Gaussian kernel and then 
transforms it for the Epanechnikov kernel using Equation (17). The bandwidths h\ and h 2 
of two kernels k x and k 2 are related as: 




(17) 



fu 
hi 



(18) 
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where li 2 (k) is the second order moment and R{k) is the roughness of the kernel k. 

[0056] The Sheather- Jones direct-plug-in method assuming a Gaussian kernel is 
described as follows. First, calculate the estimate of scale (A) for the data. It is 

preferred to use the smaller of the sample standard deviation (cr) and the normalized 
interquartile range (IQR) as a scale estimate: 



A = min[ <r,— 1 
I 1.349 ) 



[0057] Second, calculate the normal-scale approximation, 

ro= _105_ 
8 32^2 9 

where n is the sample size. Third, calculate the first approximation of the bandwidth 



(19) 



(20) 



-2/c (6) (0) 



30 



1/9 



ju 2 (K)y> s NS n 

where k: (6) (0) is the sixth derivative of the kernel function k evaluated at zero, and 
ju 2 (k) is the second-order moment for the kernel k . Here k = Gaussian kernel is 
assumed. Fourth, calculate the second approximation of the bandwidth using the fourth 
derivative of the kernel and the kernel estimator y/ 6 (g\) : 



(21) 



-2*r w (0) 
M^te 1 )" L2*i-15n_ 



-6n 



1/7 



(22) 



where: 



si = 



n~\ n 



Z Z e-" j/2 [([tu-l5] ti ,j + 45)tu-\5] 



, and 



ti,j = 



Xi — Xj 

K 8* , 



(23) 
(24) 



[0058] Finally, calculate the bandwidth using the kernel roughness function 
R(k) and the kernel estimator y/ A {g2) : 
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h = 



R(K) 



-il/5 



n 



72(2^ + 3/i) 



£2 



(25) 



where: 



Si = 



I I'"" J "([' , «-«]''^') 

_ 1=1 7=i+l 



, and 



(26) 



(27) 



v # 2 ; 

[0059] A comparison of the kernel density estimates using the first generation 
method and Sheather- Jones direct-plug-in bandwidths is presented in FIGURE 7 for a 



multi-modal density p(x) = — N 

20 



6 
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kernel density estimate using the first generation bandwidth method completely misses 
the middle peak at x = 0, while the density estimate using the SJ-dpi bandwidth is able to 
reveal the tri-modal structure of the data. Thus, FIGURE 7 shows the improvement of 
the SJ-dpi method over the first generation bandwidth calculation method. 

[0060] The difference between the actual density and the calculated kernel 
density is called the bias of the kernel density. Kernel density estimates can have large 
bias near boundaries when the data have lower or upper bounds. For the present 
embodiment on chiller vibration data, the vibration amplitudes are always positive, that 
is, the data have a lower boundary at zero. To correct the error in estimating density near 
the lower boundary (near zero), a different set of kernels may be used near the boundary. 
For kernels with a [-1, 1] support, such as the Epanechnikov kernel, boundary kernels 
may be calculated as: 



K L {u,a) 



k(u)\<_ 



{-K«<fl} 



(28) 
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where k l (•) is the boundary kernel, k (•) is a kernel function defined in Table 1 with a 
[-1, 1] support, u = (x-x i )/h anda=x/h. The coefficients v,(a)(/ = 0,1,2) are defined 
as: 



v,(a) = \u'tc(u)du (29) 



For the Epanechnikov kernel, the coefficients v/ (a) are: 

v 0 (a)=4(-a 3 +3a + 2) 
4 

v 1 ( a ) = -^(-a 4 +2fl 2 -l) (30) 
lo 

v 2 (a) = — (~3a 5 +5a 3 +2) 
2 20 

The boundary kernels are used in Equation (13) when a < 1 and the standard kernel is 
used for a >1 . The Epanechnikov boundary kernels for different values of the parameter 
a = x/h are presented in FIGURE 8. 

[0061] Because the boundary kernels have negative values, the density estimate 
using these kernels can also be negative. Thus, in order to obtain a positive density 
estimate, the negative density may be truncated to zero and then the density estimate may 
be re-scaled as: 

/**)<- r P(x) (3D 

I p(x)dx 

such that p(x) now satisfies p(x) > 0 and J£° p{x)dx = 1 . 

[0062] In addition to boundary kernels, two other methods for reducing 
boundary bias may be utilized such as data reflection and data transformation. The data 
reflection method makes the density estimate "consistent" but still produces estimates 
with a large boundary bias. Thus, boundary kernels are needed to reduce the bias at the 
boundary to the same level as that in the interior. The data transformation method is 
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more complicated and computationally intensive than the boundary kernels approach. 
Thus, the boundary kernels is a preferred method to reduce the boundary bias in kernel 
density estimation. 

[0063] A simple example illustrating the effect of boundary bias on kernel 
density estimation is presented in FIGURE 9 for the exponential density p(x) = e x . The 
exponential density has a lower boundary at x = 0 so that: 



In this example, one thousand observations were sampled from the exponential 
distribution and the data were used to estimate the kernel density. FIGURE 9 shows that 
with no boundary bias correction, the kernel density estimate near the boundary is very 
different from the true density. The estimate is better with the reflection method, while 
the boundary kernels approach produces very accurate results. 

[0064] Once the conditional probability density for the amplitude is known, it is 
used to calculate limits as follows. The limits for the vibration amplitudes may be 

estimated by first numerically integrating the estimated conditional density p(x\F = f) in 
Equation (14) using Simpson's rule to obtain the cumulative conditional density function 
P{x\F = /) . A cutoff value, /?, is also specified for the cumulative probability function. 
Then, the amplitude limit xp is estimated by calculating the inverse of the cumulative 
probability density function P(x\F = f) at a probability value of j3: 



[0065] Typically, a vibration analyst sets three levels of amplitude limits that 
are referred to as the "alert," "alarm" and "danger" limits. Common values for the 
probability cutoffs for these limits are /5 = 0.95, 0.99 and 0.995 respectively. 




(32) 



= /) = /? 



(26) 



Table 1. Probability cutoff limits for different importance levels. 
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Importance 
Level 


0 


Alert 


Alarm 


Danger 


1 (low) 


0.9995 


0.9999 


0.99995 


2 


0.999 


0.9995 


0.9999 


3 


0.995 


0.999 


0.9995 


4 


0.99 


0.995 


0.999 


5 (high) 


0.95 


0.99 


0.995 



[0066] There are special frequencies in the vibration spectrum that are called 
"diagnostic frequencies" because they correspond to the characteristic frequencies of 
different machine components. Vibration analysts commonly analyze vibrations near 
these frequencies to detect mechanical faults. Some examples of diagnostic frequencies 
for chillers are motor, compressor, gear-mesh and blade-pass frequencies. 

[0067] In the exemplary methodology, a vibration analyst can specify 
importance levels for different frequencies in the vibration spectrum. For example, for a 
Chiller A, the motor frequencies are assigned high importance while the blade-pass 
frequencies and other parts of the spectrum have relatively lower importance levels. The 
probability cutoffs are determined from the importance level of a frequency. 

[0068] A five level importance scheme may be used where each importance 
level corresponds to a different value for the probability cutoff for each of the alert, alarm 
and danger limits. According to an exemplary embodiment, an importance level of "one" 
means low importance, while a value of "five" means high importance. The probability 
cutoff values for the five importance levels are presented in Table 1. Typically, the 
diagnostic frequencies of a machine (e.g., motor, blade-pass, gear-mesh frequencies, etc.) 
have high importance levels, while other parts of the spectrum have lower importance. 
Hence, the default importance level for the entire spectrum except the diagnostic 
frequencies is set to "one" (low importance). 

[0069] Low importance levels typically have high probability cutoff values. 
These high probability cutoffs result in higher amplitude limits for frequencies that are 
not very important. Higher amplitude limits for less important frequencies will result in a 
lower number of false alarms. Thus, an importance level scheme such as described 



-19- 



Atty. Dkt. No.: 081445-0359 



herein allows the vibration analyst to focus on the important frequencies while retaining 
the ability to detect faults when unusually high amplitudes are observed at less important 
frequencies. 

[0070] Examination of the foregoing importance level scheme reveals that the 
probability cutoff limits have a sharp discontinuity at the diagnostic frequencies. For 
example, the importance level changes from "one" to "five" at the motor frequencies and 
results in an abrupt change in the probability cutoff values. This abrupt change in the 
probability cutoffs may result in a sharp discontinuity for the amplitude limit envelope at 
the diagnostic frequencies. To have a continuous transition between the cutoff limits, a 
weighting scheme is used to smooth the probability cutoff values near the diagnostic 
frequencies: 



Pf = fidef + {fidiag - fidef) X eXp 



4N bw 



(34) 



[0071] In Equation (34), /3f is the probability cutoff value for frequency 
F = f,f3def is the default probability cutoff corresponding to the lowest importance level, 
and fidiag is the probability cutoff limit for the closest diagnostic frequency fdiag . If the 
frequency /is between two diagnostic frequencies that are closer than 12* Ntw from each 
other, then linear interpolation may be used to determine the probability cutoff value 
between the two diagnostic frequencies, otherwise Equation (34) is used. This situation 
occurs while calculating the probability cutoff values between 2xmotor and electrical 
frequencies. FIGURE 10 shows the probability cutoff levels for the 0-200 Hz spectrum 
with diagnostic frequencies of half, one, two and three times the motor frequency. The 
importance levels for these diagnostic frequencies are "five," "five," "four" and "three" 
respectively. 

[0072] According to an exemplary embodiment, vibration spectra for two 
chillers (Chiller A and B) were obtained and outlier removal and density estimation were 
performed for every frequency line using methods such as described above. In addition, 
the limits were calculated for every frequency line. 
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[0073] The results for the motor and compressor vertical positions for Chiller A 
are presented in FIGURES 1 1 and 12 respectively. These figures show that although it is 
useful to detect problems at the diagnostic frequencies of motor harmonics, sub- 
harmonics and blade-pass frequencies, other frequency ranges play a role. For example, 
a high amplitude of vibration in frequency band 1-1 .5 x motor frequency may signal a 
bearing fault that will be missed if only the diagnostic frequencies were analyzed. The 
exemplary approach allows the vibration analyst to quickly compare the entire spectrum 
with the limit envelopes to detect problems in any part of the spectrum. 

[0074] It is important to note that the above-described preferred embodiments 
are illustrative only. Although the invention has been described in conjunction with 
specific embodiments thereof, those skilled in the art will appreciate that numerous 
modifications are possible without materially departing from the novel teachings and 
advantages of the subject matter described herein. Accordingly, these and all other such 
modifications are intended to be included within the scope of the present invention as 
defined in the appended claims. The order or sequence of any process or method steps 
may be varied or re-sequenced according to alternative embodiments. In the claims, any 
means-plus-function clause is intended to cover the structures described herein as 
performing the recited function and not only structural equivalents but also equivalent 
structures. Other substitutions, modifications, changes and omissions may be made in the 
design, operating conditions and arrangements of the preferred and other exemplary 
embodiments without departing from the spirit of the present invention. 
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